#ifndef _FE_VEC_MAT_H_


#ifdef VECMAT_MPI
#include <mpi.h>
#include "fe_mpi_tags.h"
#endif

class Matrix;
class Vector;

//! CompressedSparseRow Matrix
struct CSRMat {
	CSRMat();
	~CSRMat();
	CSRMat operator= ( const CSRMat & C );
	CSRMat( const CSRMat & C );
	int	* ap, * ai, dim, nnz;
	double * a;
	bool IS_SYMMETRIC;
};

void		FreeDblArray3D( double *** &A, const int n, const int m );

void		FreeDblArray1D(	double * &x); /** free vector */
void		NewDblArray1D(	double * &x, const int n ); /** return new n vector */
double*		NewDblArray1D(	const int n ); /** return new n vector */

void		PrintFullMatrix(double ** A, const int n, const int m=-1);
double**	FullMatrix( const CSRMat & K);
void		FreeDblArray2D(	double ** &A, const int n); /** free matrix */
double**	NewDblArray2D(	const int n, const int m ); /** return new n x m matrix */
double**	MatCopy(	double ** A, const int n, const int m ); /** return copy of n x m matrix */
double**	IdentityMatrix(	const int n ); /** return new n x n identity matrix */

void		GetTransposedMat(	double **&A, const int n, const int m ); /** Transpost matrix: A -> A^T */
double**	TransposeMat(		double **A, const int n, const int m ); /** return transposed matrix */

void		ScalVecMul(		const double a, double * b, const int n); /** Vector times scalar b -> a*b */
void		AddVectors(		double * z, const double * x, const double * y, const int n, const double a = 1.); /** sum of two vectors: z = x
+y */
double*		AddVectors(		const double * x, const double * y, const int n, const double a = 1.);  /** returns sum of two vectors */

double*		VecLinComb( double ** X, double * a, const int n, const int m ); /** return Vector linear combination */
void		EraseVector( double * x, const int N, const double z = 0. );
void		EraseMatrix( double ** x, const int N, const int M, const double z = 0. );

double VectorNorm(const double * x, const int n = 3, const int p = 2);
double DeviNorm(const double * x );

double		VecVecMul( 		const double * x, const double * y, const int n ); /** Vector Vector scalar multiplication */
double*		MatVecMul(		double ** A, const double * x, const int n, const int m); /** return Vector Matrix Product Ax */
double*		MatTVecMul(		double ** A, const double * x, const int n, const int m); /** return Vector Matrix Product Ax */
void		MatVecMul(		double ** A, const double * x, double * y, const int n, const int m); /** return Vector Matrix Product Ax */
void		MatTVecMul(		double ** A, const double * x, double * y, const int n, const int m); /** return Vector Matrix Product Ax */
void		MatAddMul(		double ** &A, double ** B,	const int n, const int m, const double z = 1. );
void		MatMultiply(double ** &A, const int n, const int m, const double z = 1.);
void 		MatMatMul(		double ** A, double ** B,	double ** C, const int n, const int m, const int o, const int p);
double**	MatMatMul(		double ** A, double ** B,	const int n, const int m, const int o, const int p);
	/** return Matrix Matrix Product A*B, A in R ^ n x m, B in R ^ m x o */
void		MatTMatMul(		double ** A, double ** B, double ** C,	const int n, const int m, const int o, const int p);
double**	MatTMatMul(		double ** A, double ** B,	const int n, const int m, const int o, const int p);
	/** return Matrix Matrix Product A^T*B, A in R ^ m x n, B in R ^ m x o */
void		MatMatTMul(		double ** A, double ** B, double ** C,	const int n, const int m, const int o, const int p);
double**	MatMatTMul(		double ** A, double ** B,	const int n, const int m, const int o, const int p);
	/** return Matrix Matrix Product A*B^T, A in R ^ n x m, B in R ^ o x m */
void		MatTMatTMul(		double ** A, double ** B, double ** C,	const int n, const int m, const int o, const int p);
double**	MatTMatTMul(		double ** A, double ** B,	const int n, const int m, const int o, const int p);
	/** return Matrix Matrix Product A^T*B^T, A in R ^ m x n, B in R ^ o x m */

void		SolveLinearSystems( double ** A, double ** y, double ** x, const int N, const int M);
void		SolveLinearSystem(double ** A, const double * y, double * x, const int N);
void		InvertMatrix(	double ** A, double ** & invA, const int N ); /** comput inverse of A */
double**	InvertMatrix(	double ** A, const int N ); /** return A^-1 */

double		Det3x3(			double ** A ); /** return det(A) for A in R^3x3 */

/* ************************************************************************ */
/* ************************************************************************ */
/* ************************************************************************ */

Vector		NewVector( const int n );

void		FreeVector( Vector & x );
void		FreeMatrix( Matrix & M );
Matrix		NewMatrix( const int n, const int m = 0);
Matrix		MatCopy(	const Matrix & B );
void		MatCopy(	Matrix & A, const Matrix & B );
Matrix		MIdentityMatrix( const int n = 3 );

void		GetTransposedMat(	Matrix & A );
Matrix		TransposeMat(		const Matrix & A );

void		ScalVecMul(		const double a, Vector & b );
void		AddVectors(		Vector & x, const Vector & y, const Vector & z, const double a=1. );
Vector		AddVectors(		const Vector & x, const Vector & y, const double a=1. );
void		EraseVector( Vector & x, const double z = 0. );
void		EraseMatrix( Matrix & M, const double z = 0. );

double		VectorNorm(const Vector & x, const int p = 2);

double		VecVecMul( 		const Vector & x, const Vector & y ); /** Vector Vector scalar multiplication */
Vector		MatVecMul(		const Matrix & A, const Vector & x ); /** return Vector Matrix Product Ax */
Vector		MatTVecMul(		const Matrix & A, const Vector & x ); /** return Vector Matrix Product Ax */
void		MatVecMul(		Vector & y, const Matrix & A, const Vector & x ); /** return Vector Matrix Product Ax */
void		MatTVecMul(		Vector & y, const Matrix & A, const Vector & x );/** return Vector Matrix Product Ax */
void		MatAddMul(		Matrix & A, const Matrix & B, const double z = 1. );
void		MatMultiply(	Matrix & A, const double z = 1.);
void 		MatMatMul(		const Matrix & A, const Matrix & B,	Matrix & C );
Matrix		MatMatMul(		const Matrix & A, const Matrix & B );
	/** return Matrix Matrix Product A*B, A in R ^ n x m, B in R ^ m x o */
void		MatTMatMul(		const Matrix & A, const Matrix & B, Matrix & C );
Matrix		MatTMatMul(		const Matrix & A, const Matrix & B );
	/** return Matrix Matrix Product A^T*B, A in R ^ m x n, B in R ^ m x o */
void		MatMatTMul(		const Matrix & A, const Matrix & B, Matrix & C );
Matrix		MatMatTMul(		const Matrix & A, const Matrix & B );
	/** return Matrix Matrix Product A*B^T, A in R ^ n x m, B in R ^ o x m */
void		MatTMatTMul(		const Matrix & A, const Matrix & B, Matrix & C );
Matrix		MatTMatTMul(		const Matrix & A, const Matrix & B );
	/** return Matrix Matrix Product A^T*B^T, A in R ^ m x n, B in R ^ o x m */

void		SolveLinearSystems( const Matrix & A, const Matrix & Y, Matrix & X);
void		SolveLinearSystem( const Matrix & A, const Vector & y, Vector & x);
void		InvertMatrix(	const Matrix & A, Matrix & invA ); /** comput inverse of A */
Matrix		InvertMatrix(	const Matrix & A ); /** return A^-1 */

double		Det3x3(			Matrix & A ); /** return det(A) for A in R^3x3 */


/* ************************************************************************ */
/* ************************************************************************ */
/* ************************************************************************ */

//! Simple datatype for a dense matrix (always non-symmetric)
class Matrix {
// 	Matrix() : A(0),m(0),n(0) { };
// 	protected:
		//! Array containing the double precision values of the matrix; Size m x n
		double	**A; /* Attention: PSEUDO 2D ARRAY */
		bool	is_diag, is_dummy;
	public:
		//! Number of rows
		int		m;
		//! Number of columns
		int		n;
		Matrix();
		Matrix(const int N, const bool del = true );
		Matrix(const int M, const int N, const bool del = true );
		Matrix(const Matrix & B);
		Matrix(double * x, const int m, const int n = -1);
		~Matrix();
		void print() const;
		void Identity(); /* Set Matrix to identity matrix */
		inline bool IsDiagonal() const { return is_diag; }
		inline void SetDiagonal( const bool diag ) { is_diag = diag; }

#ifdef VECMAT_MPI
		void MPIBroadcast(const int root = 0);
		void MPISendToFrom( const int to, const int from = 0);
#endif
		double *	GetArray1D() const {	return A[0];	}
		double **	GetArray2D() const {	return A;		}
		//! Dense vector matrix multiplication
		inline double& operator() (const int i, const int j) { return A[0][i*n+j]; }; // return A[i][j]; };
		inline const double& operator() (const int i, const int j) const { return A[0][i*n+j]; }; // return A[i][j]; };
		void resize(const int newM, const int newN, const bool del = false );
		void ReadDataFromFile( FILE * F, const int M, const int N );
		inline double * operator() (const int i) const { return A[i]; }
		inline double * operator[] (const int i) const { return A[i]; }
// 		inline const double * operator() const (const int i) { return A+i*n; }
		Matrix& operator=(const Matrix & B);
		Matrix& operator=(const double a);
		Matrix& operator*=(const double a);
		Matrix& operator/=(const double a);
		Matrix& operator+=(const Matrix & B);
		Matrix& operator-=(const Matrix & B);
	// 	Vector operator*(const Vector & y) const;
		//! Sum of dense matrices
		Matrix operator+(const Matrix & B) const;
		//! Difference of dense matrices
		Matrix operator-(const Matrix & B) const;
		void SymEigenValVec( Matrix & EV, Vector & lambda );
		void transpose();
};

/* ************************************************************************ */
/* ************************************************************************ */
/* ************************************************************************ */

//! Simple datatype for a dense vector
class Vector {
	public:
		Vector();
		Vector(const Vector & y);
		Vector(const int N, const bool del = true );
		~Vector();
#ifdef VECMAT_MPI
		void MPIBroadcast(const int root = 0);
#endif
		//! Array containing the double precision values of the vector
		double	*X;
		//! dimension of the vector
		int		n;
		void print() const;

		void resize(const int N, const bool del = false );
		Vector& operator= (const double a);
		Vector& operator= (const Vector & y);
		Vector& operator+= (const Vector & y);
		Vector& operator-= (const Vector & y);
		Vector& operator*= (const double a);
		void ReadDataFromFile( FILE * F, const int N );
		//! inner product of two vectors
		inline double& operator() (const int i) { return X[i]; };
		inline const double& operator() (const int i) const { return X[i]; };
		inline double& operator[] (const int i) { return X[i]; };
		inline const double& operator[] (const int i) const { return X[i]; };
		double operator*(const Vector & y) const;
		//! sum of vectors
		friend Vector operator-(const Vector & x);
		Vector operator+(const Vector & y) const;
		//! difference of vectors
		Vector operator-(const Vector & y) const;
};

void BLAS_MatMatMul(	const Matrix & A, const Matrix & B, Matrix & C );
void BLAS_MatTMatMul(	const Matrix & A, const Matrix & B, Matrix & C );
void BLAS_MatMatTMul(	const Matrix & A, const Matrix & B, Matrix & C );
void BLAS_MatTMatTMul(	const Matrix & A, const Matrix & B, Matrix & C );
/* ************************************************************************ */
/* ************************************************************************ */
/* ************************************************************************ */

//! multiply dense matrix by scalar value
Vector operator*(const Vector & y, const double a) ;
Vector operator*(const double a,   const Vector & y) ;
Vector operator/(const Vector & y, const double a) ;
Matrix operator*(const Matrix & A, const double a) ;
Matrix operator/(const Matrix & A, const double a) ;
Matrix operator*(const Matrix & A, const Matrix & B) ;
Matrix operator*(const double a, const Matrix & A) ;
//! dense vector matrix multiplication
Vector operator*(const Matrix & A, const Vector & y) ;
Vector operator*(const Vector & y, const Matrix & A) ;
Vector operator/(const Matrix & A, const Vector & y) ;

/* ************************************************************************ */
/* ************************************************************************ */
/* ************************************************************************ */

#define _FE_VEC_MAT_H_


#endif /* _FE_VEC_MAT_H_ */


